library(TTR)
library(quantmod)

library(ggplot2); library(readr); library(dplyr); library(tidyverse);
library(gridExtra);library(ggExtra); library(knitr);
library(caret); library(gridExtra); library(lubridate); library(stringr);
library(PerformanceAnalytics); library(xts); library(tidyquant); library(tidyverse); library(data.table); library(chron)

Mass shootings data

us_mass <- read.csv("US_Mass_Shootings_May_24_2022.csv")

#Select those variables that are needed for our analysis
us_mass02 <- us_mass %>% select(location, date, fatalities, injured, total_victims)
str(us_mass02)
## 'data.frame':    128 obs. of  5 variables:
##  $ location     : chr  "Uvalde, Texas" "Buffalo, New York" "Sacramento, California" "Oxford, Michigan" ...
##  $ date         : chr  "5/24/22" "5/14/22" "2/28/22" "11/30/21" ...
##  $ fatalities   : int  15 10 4 4 9 8 4 10 8 4 ...
##  $ injured      : chr  "-" "3" "0" "7" ...
##  $ total_victims: chr  "-" "13" "4" "11" ...
# change date format to standardize dates
us_mass02 <- us_mass02 %>% mutate(date = mdy(date)) %>%
                #Shooting was/is in a holiday date is move to next date
                mutate(date = ifelse(is.holiday(date), date+1,date)) %>%
                # Moving any dates to the following Monday
                mutate(weekday=wday(date)) %>%
                mutate(date1 = ifelse(weekday==1,(date+1), ifelse(weekday==7, (date+2), (date)))) %>%
                mutate(date = ifelse(is.holiday(date), date+1,date)) %>%
                mutate(date = as.Date(date1)) %>%
                mutate(year = year(date)) %>%
               
                mutate(state = ifelse(location=="Washington, D.C.", "Washington, D.C.",str_extract(location,
                                                "[\\s]*[A-Z][a-z]+[\\s]*[[A-Z][a-z]+]*$"))) %>%
                mutate(fatalities = as.integer(fatalities)) %>%
                mutate(injured = as.integer(injured)) %>%
                mutate(total_victims = as.integer(total_victims)) %>%
                select(location, state, date, year, fatalities, injured, total_victims)
## Warning in mask$eval_all_mutate(quo): NAs introduced by coercion

## Warning in mask$eval_all_mutate(quo): NAs introduced by coercion
## Update Uvalde, Texas data
us_mass02$fatalities[us_mass02$location=="Uvalde, Texas"] <- as.integer(22)
us_mass02$injured[us_mass02$location=="Uvalde, Texas"] <- as.integer(17)
us_mass02$total_victims[us_mass02$location=="Uvalde, Texas"] <- as.integer(39)

## Add Highland Park, Illinois mass shooting  
us_mass02 <- us_mass02 %>% add_row(location = "Highland Park, Illinois", state = "Illinois", date = mdy(07042022), year= as.integer(2022), fatalities = as.integer(8), injured=as.integer(29), total_victims=as.integer(37), .before=1)




summary(us_mass02)
##    location            state                date                 year     
##  Length:129         Length:129         Min.   :1982-08-20   Min.   :1982  
##  Class :character   Class :character   1st Qu.:2001-02-05   1st Qu.:2001  
##  Mode  :character   Mode  :character   Median :2013-06-07   Median :2013  
##                                        Mean   :2009-10-08   Mean   :2009  
##                                        3rd Qu.:2018-02-14   3rd Qu.:2018  
##                                        Max.   :2022-07-04   Max.   :2022  
##    fatalities        injured       total_victims   
##  Min.   : 3.000   Min.   :  0.00   Min.   :  3.00  
##  1st Qu.: 4.000   1st Qu.:  1.00   1st Qu.:  7.00  
##  Median : 6.000   Median :  3.00   Median : 10.00  
##  Mean   : 8.093   Mean   : 11.69   Mean   : 19.78  
##  3rd Qu.: 9.000   3rd Qu.: 10.00   3rd Qu.: 18.00  
##  Max.   :58.000   Max.   :546.00   Max.   :604.00

Group mass shooting data by year and state

Group mass shooting data by state and year

usmass_yearly <- us_mass02 %>% group_by(year) %>%
        summarise( fatalities=sum(fatalities, na.rm = TRUE), injured=sum(injured, na.rm = TRUE), total_victims=sum(total_victims, na.rm = TRUE)) 

shootings_yrl <- us_mass02 %>% count(year)
colnames(shootings_yrl) <- c("year", "mass_shootings")
usmass_yearly <- merge(usmass_yearly, shootings_yrl, by="year")
usmass_yearly$intensity <- usmass_yearly$total_victims/usmass_yearly$mass_shootings
tail(usmass_yearly)
##    year fatalities injured total_victims mass_shootings intensity
## 33 2017        117     587           704             11 64.000000
## 34 2018         80      70           150             12 12.500000
## 35 2019         73     112           185             10 18.500000
## 36 2020          9       0             9              2  4.500000
## 37 2021         43      16            59              6  9.833333
## 38 2022         44      49            93              4 23.250000
usmass_state <- us_mass02 %>% group_by(state) %>%
        summarise( fatalities=sum(fatalities, na.rm = TRUE), injured=sum(injured, na.rm = TRUE), total_victims=sum(total_victims, na.rm = TRUE))

shootings_state <- us_mass02 %>% count(state)
colnames(shootings_state) <- c("state", "mass_shootings")
usmass_state <- merge(usmass_state, shootings_state, by="state")
usmass_state <- usmass_state %>% arrange(desc(mass_shootings))
head(usmass_state)
##         state fatalities injured total_victims mass_shootings
## 1  California        157     160           317             23
## 2     Florida        126     109           235             12
## 3       Texas        152     183           335             12
## 4    Colorado         48     104           152              7
## 5  Washington         37      28            65              7
## 6    New York         40      28            68              5

Exploratory plotting

US Mass Shooting total victims per State

Plot of total fatalities grouped by state

subset_usmass_state1 <- usmass_state %>% arrange(desc(total_victims)) %>% slice(1:12)
temp_masssate <- subset_usmass_state1[,c("state", "injured", "fatalities")] %>% reshape2::melt()
## Using state as id variables
usshoot_state_plot <- ggplot(data = temp_masssate, aes(y=reorder(state, (value)), x=value, fill=variable))
usshoot_state_plot <- usshoot_state_plot + geom_col(position="stack", stat="identity")
## Warning: Ignoring unknown parameters: stat
usshoot_state_plot <- usshoot_state_plot + scale_x_continuous(breaks = seq(0, 700, by = 100))
usshoot_state_plot <- usshoot_state_plot + labs(title="US Mass Shooting Total Victims per State",
                                                subtitle = "Kaggle US Mass Shootings",
                                                caption = "Team 89")
usshoot_state_plot <- usshoot_state_plot + ylab("State") + xlab("Total Victims")
usshoot_state_plot <- usshoot_state_plot + theme(plot.title = element_text(size = 15))
usshoot_state_plot <- usshoot_state_plot + theme(axis.title = element_text(size = 15)) 
usshoot_state_plot <- usshoot_state_plot + theme(axis.text = element_text(size = 15))
usshoot_state_plot <- usshoot_state_plot + theme(legend.title = element_blank())
usshoot_state_plot <- usshoot_state_plot + theme(legend.position = "right")




usshoot_state_plot

ggsave("89plot01.png")
## Saving 7 x 5 in image

US Mass Shooting total shootings per State

Plot of number of mass shootings grouped by state

#temp_masssate <- usmass_state[,c("state", "injured", "fatalities")] %>% reshape2::melt()

subset_usmass_state2 <- usmass_state %>% arrange(desc(mass_shootings)) %>% slice(1:12)
        
usshoot_state_plot1 <- ggplot(data = subset_usmass_state2, aes(y=reorder(state, mass_shootings), x=mass_shootings))
usshoot_state_plot1 <- usshoot_state_plot1 + geom_col(position="stack", stat="identity", fill="red")
## Warning: Ignoring unknown parameters: stat
usshoot_state_plot1 <- usshoot_state_plot1 + scale_x_continuous(breaks = seq(0, 30, by = 5))
usshoot_state_plot1 <- usshoot_state_plot1 + labs(title="US Total of Mass Shootings per State",
                                                subtitle = "Kaggle US Mass Shootings",
                                                caption = "Team 89")
usshoot_state_plot1 <- usshoot_state_plot1 + ylab("State") + xlab("Total Mass Shootings")
usshoot_state_plot1 <- usshoot_state_plot1 + theme(plot.title = element_text(size = 15))
usshoot_state_plot1 <- usshoot_state_plot1 + theme(axis.title = element_text(size = 15)) 
usshoot_state_plot1 <- usshoot_state_plot1 + theme(axis.text = element_text(size = 15))
usshoot_state_plot1 <- usshoot_state_plot1 + theme(legend.title = element_blank())
usshoot_state_plot1 <- usshoot_state_plot1 + theme(legend.position = "right")

usshoot_state_plot1

ggsave("89plot02.png")
## Saving 7 x 5 in image

Annual Mass Shootings

Plot of number of mass shooting per year

usshoot_plot1 <- ggplot()+ geom_point(data = usmass_yearly, aes(x=year, y=mass_shootings), color="red")
usshoot_plot1 <- usshoot_plot1 + geom_smooth(data = usmass_yearly, aes(x=year, y=mass_shootings), method="lm")
usshoot_plot1 <- usshoot_plot1 + labs(title = "Annual Mass Shootings", 
                                      subtitle = "Kaggle US Mass Shootings", 
                                      caption = "Team 89")

usshoot_plot1
## `geom_smooth()` using formula 'y ~ x'

ggsave("89plot03.png")
## Saving 7 x 5 in image
## `geom_smooth()` using formula 'y ~ x'

Annual Fatalities

Plot of annual fatalities due to mass shootings

usshoot_plot1 <- ggplot()+ geom_point(data = usmass_yearly, aes(x=year, y=fatalities), color="red")
usshoot_plot1 <- usshoot_plot1 + geom_smooth(data = usmass_yearly, aes(x=year, y=fatalities), method="lm")
usshoot_plot1 <- usshoot_plot1 + labs(title = "Annual Mass Shooting Fatalities", 
                                      subtitle = "Kaggle US Mass Shootings", 
                                      caption = "Team 89")

usshoot_plot1
## `geom_smooth()` using formula 'y ~ x'

ggsave("89plot04.png")
## Saving 7 x 5 in image
## `geom_smooth()` using formula 'y ~ x'

Annual Fatality Intesity

Number of victims per shooting on function of year

usshoot_plot1 <- ggplot()+ geom_point(data = usmass_yearly, aes(x=year, y=intensity), color="red")
usshoot_plot1 <- usshoot_plot1 + geom_smooth(data = usmass_yearly, aes(x=year, y=intensity), method="lm")
usshoot_plot1 <- usshoot_plot1 + labs(title = "Number of Victims per mass shooting", 
                                      subtitle = "Kaggle US Mass Shootings", 
                                      caption = "Team 89")
usshoot_plot1 <- usshoot_plot1 + ylab("Victims per mass shooting")
usshoot_plot1
## `geom_smooth()` using formula 'y ~ x'

ggsave("89plot05.png")
## Saving 7 x 5 in image
## `geom_smooth()` using formula 'y ~ x'

Obtain historical stock price of US gun manufactures.

Obtain stock history for the companies that are been evaluated using tq_get function

date_str <- "2009-10-01"
#stocks <- c("swbi")
stocks <- c("rgr", "swbi")
#stocks <- c("rgr", "swbi", "oln", "npk", "bgfv", "axon" )

gun_co_stocks <- stocks %>% tq_get(get="stock.prices", from=date_str, to=Sys.Date())
head(gun_co_stocks)
## # A tibble: 6 × 8
##   symbol date        open  high   low close volume adjusted
##   <chr>  <date>     <dbl> <dbl> <dbl> <dbl>  <dbl>    <dbl>
## 1 rgr    2009-10-01  12.9  13.0  12.4  12.5 252500     8.38
## 2 rgr    2009-10-02  12.4  12.7  11.8  12.3 254000     8.25
## 3 rgr    2009-10-05  12.3  12.5  12.2  12.2 120300     8.20
## 4 rgr    2009-10-06  12.2  12.8  12.2  12.5 111300     8.41
## 5 rgr    2009-10-07  12.5  12.8  12.4  12.6  99300     8.45
## 6 rgr    2009-10-08  12.7  13.0  12.5  12.8 144400     8.60

Plotting adjusted stock price using ggplot2

The historical stock price for analyzed public companies, since October 2009

gun_stocks_plot <- ggplot(data = gun_co_stocks, aes(x=date, y=adjusted, group=symbol))
gun_stocks_plot <- gun_stocks_plot + geom_line(aes(color=symbol))
gun_stocks_plot <- gun_stocks_plot + theme(legend.position = "bottom") 
gun_stocks_plot <- gun_stocks_plot + labs(title = "Adjusted Stock Price US Gun Manufactures", 
                                          subtitle = "Source: Yahoo Finance", 
                                          caption = "team 89")
gun_stocks_plot <- gun_stocks_plot + ylab("Adjusted Stock Price")
gun_stocks_plot <- gun_stocks_plot + scale_x_date(date_breaks = "2 years", date_labels = "%Y")

gun_stocks_plot

ggsave("89plot06B.png")
## Saving 7 x 5 in image

Filter mass shootings that happened after 2010-01-01, then Obtain mass shooting date vector

Generate a vector with dates of mass shootings so we can use it in the following exploratory plots

us_mass03 <- us_mass02 %>% filter(date>="2010-01-01" & total_victims>=7)
us_mass_dates <- us_mass03[,3]

Plotting adjusted stock price and mass shooting dates using ggplot2

Generate plots including adjusted stock price and dates of mass shootings (vertical dash lines)

gun_stocks_plot <- ggplot(data = gun_co_stocks, aes(x=date, y=adjusted, group=symbol))
gun_stocks_plot <- gun_stocks_plot + geom_line(aes(color=symbol), size=0.75)
gun_stocks_plot <- gun_stocks_plot + theme(legend.position = "bottom") 
gun_stocks_plot <- gun_stocks_plot + theme(legend.position = "bottom", legend.key.width=unit(0.1,"npc"))

## Vertical lines to show mass shootings
gun_stocks_plot <- gun_stocks_plot + geom_vline(data=us_mass03, aes(xintercept=date), linetype=3)

gun_stocks_plot <- gun_stocks_plot + labs(title = "Adjusted Stock Price US Gun Manufactures", 
                                          subtitle = "Mass Shootings with seven or more total victims", 
                                          caption = "Source: Yahoo finance and Kaggle's Mass Shooting in United States (2018-2022)")
gun_stocks_plot <- gun_stocks_plot + ylab("Adjusted Stock Price")
gun_stocks_plot <- gun_stocks_plot + scale_x_date(date_breaks = "2 years", date_labels = "%Y")

gun_stocks_plot

ggsave("89plot07B.png")
## Saving 7 x 5 in image

Get stock returns

Obtain the returns for each of companies using tq_transmute function from tidyquant package

gun_co_returns <- gun_co_stocks %>% group_by(symbol) %>% tq_transmute(select = adjusted, 
                                                                      mutate_fun = periodReturn, 
                                                                      period="daily", 
                                                                      col_rename = "Ra")

head(gun_co_returns)
## # A tibble: 6 × 3
## # Groups:   symbol [1]
##   symbol date             Ra
##   <chr>  <date>        <dbl>
## 1 rgr    2009-10-01  0      
## 2 rgr    2009-10-02 -0.0160 
## 3 rgr    2009-10-05 -0.00570
## 4 rgr    2009-10-06  0.0254 
## 5 rgr    2009-10-07  0.00559
## 6 rgr    2009-10-08  0.0175

Plot of gun company returns

Generate a plot showing return of each of the companies and dates of mass shootings (vertical dash lines)

gun_returns_plot <- ggplot(data = gun_co_returns, aes(x=date, y=Ra, group=symbol))
gun_returns_plot <- gun_returns_plot + geom_line(aes(color=symbol))
gun_returns_plot <- gun_returns_plot + theme(legend.position = "bottom") 
gun_returns_plot <- gun_returns_plot + theme(legend.position = "bottom", legend.key.width=unit(0.1,"npc"))

## Vertical lines to show mass shootings
gun_returns_plot <- gun_returns_plot + geom_vline(data=us_mass03, aes(xintercept=date), linetype=3)

gun_returns_plot <- gun_returns_plot + labs(title = "Returns of US Gun Manufactures", 
                                            subtitle = "Mass Shootings with seven or more total victims", 
                                            caption = "Source: Yahoo finance and Kaggle's Mass Shooting in United States (2018-2022)")

gun_returns_plot <- gun_returns_plot + ylab("Ra") + ylim(-0.2, 0.2)
gun_returns_plot <- gun_returns_plot + scale_x_date(date_breaks = "2 years", date_labels = "%Y")

gun_returns_plot

ggsave("89plot08B.png")
## Saving 7 x 5 in image

Obtain market value (capitalization) for each company and create market value weighted index weight

rvest https://cran.r-project.org/web/packages/rvest/rvest.pdf package is used to harvest (scrap) summary page of each company in Yahoo Finance. Once the market values are obtain a weights are determine to set up portafolio

library(rvest)
## Warning: package 'rvest' was built under R version 4.2.1
## 
## Attaching package: 'rvest'
## The following object is masked from 'package:readr':
## 
##     guess_encoding
# Define stock name
#stocks <- c("rgr", "swbi")
#stocks <- c("rgr", "swbi", "oln", "npk", "bgfv", "axon" )

Mkt_Caps <- c()
for(s in stocks){
        # Extract and transform data
        temp_Co_Metrics <- paste0("https://finance.yahoo.com/quote/", s) %>% 
                           read_html() %>%  html_table() %>% map_df(bind_cols) %>% 
                           t() %>% as_tibble()

        temp_Mkt_Cap <- temp_Co_Metrics[2,9]
        temp_Cap_Unit <- str_extract(temp_Mkt_Cap, "[A-Z]")
        temp_Cap_Amount <- as.numeric(str_extract(temp_Mkt_Cap, "\\d*\\.\\d*"))
        temp_Cap_Value <- case_when(temp_Cap_Unit== "T" ~ temp_Cap_Amount*1000*1000,
                                    temp_Cap_Unit== "B" ~ temp_Cap_Amount*1000,
                                    temp_Cap_Unit=="M" ~ temp_Cap_Amount*1,
                                    is.na(temp_Cap_Unit) ~ temp_Cap_Amount/(1000*1000))
        Mkt_Caps <- append(Mkt_Caps, temp_Cap_Value)
        }
## Warning: The `x` argument of `as_tibble.matrix()` must have unique column names if `.name_repair` is omitted as of tibble 2.0.0.
## Using compatibility `.name_repair`.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.
gun_wts <- Mkt_Caps/sum(Mkt_Caps)
#gun_wts <- c(1/6, 1/6, 1/6, 1/6, 1/6, 1/6)

co_summary <- data.frame(stocks, Mkt_Caps, gun_wts)
co_summary              
##   stocks Mkt_Caps   gun_wts
## 1    rgr 1112.000 0.6413195
## 2   swbi  621.925 0.3586805

Create a portfolio of gun manufactures

A data set for the portafolio was created as follows:

  • Portfolio return, Ra, was calculated using weights determined previously and using tq_portfolio function
  • Simple moving average using tq_mutate function and mutate_fun “SMA”
  • Return drawdown using Drawdowns function. Then for each observation a new variables were created for the future drawdown at different days after using lead function
  • Portfolio rolling average using frollmean
  • Portfolio return after 2, 3, 5, 10, 15, and 20 days. Momentum is defined as \(Ra_t-Ra_0\) using lead function
  • Dataset was filter for dates after January 1, 2022
gun_co_portfolio_returns <- gun_co_returns %>% 
                                tq_portfolio(assets_col = symbol, 
                                             returns_col = Ra, 
                                             weights = gun_wts, 
                                             col_rename = "Ra")
## Warning: `spread_()` was deprecated in tidyr 1.2.0.
## Please use `spread()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.
gun_co_portfolio_returns <- gun_co_portfolio_returns %>% tq_mutate(select = Ra, mutate_fun = SMA, n=2) %>%
                                                        rename(SMA.2day = SMA) %>%
                                                        tq_mutate(select = Ra, mutate_fun = SMA, n=15) %>%
                                                        rename(SMA.15day = SMA)
                                                        


gun_co_portfolio_returns <- gun_co_portfolio_returns %>% mutate(gunDD = Drawdowns(Ra)) %>%
                                                        mutate(gunDD=as.numeric(gunDD)) %>%
                                                        mutate(gunDD_2day=lead(gunDD,2)) %>%
                                                        mutate(gunDD_3day=lead(gunDD,3)) %>%
                                                        mutate(gunDD_5day=lead(gunDD,5)) %>%
                                                        mutate(gunDD_10day=lead(gunDD,10)) %>%
                                                        mutate(gunDD_15day=lead(gunDD,15)) %>%
                                                        mutate(gunDD_20day=lead(gunDD,20))
## Warning in merge.zoo(fx, .xts(, .index(x))): Index vectors are of different
## classes: integer POSIXct
gun_co_portfolio_returns <- gun_co_portfolio_returns %>% mutate(Ra_2day=frollmean(Ra,2)) %>%
                                                        mutate(Ra_3day=frollmean(Ra,3)) %>% 
                                                        mutate(Ra_5day=frollmean(Ra,5)) %>%
                                                        mutate(Ra_10day=frollmean(Ra,10))


gun_co_portfolio_returns <- gun_co_portfolio_returns %>% mutate(Mom_2day=lead(Ra,2) - Ra) %>%
                                                        mutate(Mom_3day=lead(Ra,3) - Ra) %>%
                                                        mutate(Mom_5day=lead(Ra,5) - Ra) %>%
                                                        mutate(Mom_10day=lead(Ra,10) - Ra) %>%
                                                        mutate(Mom_15day=lead(Ra,15) - Ra) %>%
                                                        mutate(Mom_20day=lead(Ra,20) - Ra)

gun_co_portfolio_returns <- gun_co_portfolio_returns %>% filter(date>="2010-01-01")

head(gun_co_portfolio_returns,20)
## # A tibble: 20 × 21
##    date              Ra  SMA.2day SMA.15day  gunDD gunDD_2day gunDD_3day
##    <date>         <dbl>     <dbl>     <dbl>  <dbl>      <dbl>      <dbl>
##  1 2010-01-04  0.0398    0.00672   0.00234  -0.209     -0.174     -0.173
##  2 2010-01-05  0.0255    0.0326    0.00390  -0.189     -0.173     -0.179
##  3 2010-01-06  0.0177    0.0216    0.00502  -0.174     -0.179     -0.177
##  4 2010-01-07  0.00170   0.00969   0.00560  -0.173     -0.177     -0.187
##  5 2010-01-08 -0.00690  -0.00260   0.00413  -0.179     -0.187     -0.182
##  6 2010-01-11  0.00154  -0.00268   0.00489  -0.177     -0.182     -0.181
##  7 2010-01-12 -0.0123   -0.00539   0.00340  -0.187     -0.181     -0.157
##  8 2010-01-13  0.00697  -0.00268   0.00346  -0.182     -0.157     -0.141
##  9 2010-01-14  0.000996  0.00398   0.00385  -0.181     -0.141     -0.164
## 10 2010-01-15  0.0294    0.0152    0.00590  -0.157     -0.164     -0.178
## 11 2010-01-19  0.0186    0.0240    0.00720  -0.141     -0.178     -0.183
## 12 2010-01-20 -0.0270   -0.00421   0.00633  -0.164     -0.183     -0.180
## 13 2010-01-21 -0.0159   -0.0215    0.00484  -0.178     -0.180     -0.181
## 14 2010-01-22 -0.00617  -0.0110    0.00317  -0.183     -0.181     -0.181
## 15 2010-01-25  0.00366  -0.00125   0.00517  -0.180     -0.181     -0.193
## 16 2010-01-26 -0.00193   0.000863  0.00239  -0.181     -0.193     -0.212
## 17 2010-01-27  0.000884 -0.000523  0.000749 -0.181     -0.212     -0.197
## 18 2010-01-28 -0.0145   -0.00683  -0.00140  -0.193     -0.197     -0.192
## 19 2010-01-29 -0.0242   -0.0194   -0.00312  -0.212     -0.192     -0.184
## 20 2010-02-01  0.0194   -0.00237  -0.00137  -0.197     -0.184     -0.207
## # … with 14 more variables: gunDD_5day <dbl>, gunDD_10day <dbl>,
## #   gunDD_15day <dbl>, gunDD_20day <dbl>, Ra_2day <dbl>, Ra_3day <dbl>,
## #   Ra_5day <dbl>, Ra_10day <dbl>, Mom_2day <dbl>, Mom_3day <dbl>,
## #   Mom_5day <dbl>, Mom_10day <dbl>, Mom_15day <dbl>, Mom_20day <dbl>

Gun Portfolio statistics

Gun portfolio statistics using table.Stats from PerformaceAnalytics package https://www.rdocumentation.org/packages/PerformanceAnalytics/versions/2.0.4/topics/table.Stats

gun_portfolio_stats <- table.Stats(gun_co_portfolio_returns$Ra)
gun_portfolio_stats
##                          
## Observations    3160.0000
## NAs                0.0000
## Minimum           -0.1771
## Quartile 1        -0.0107
## Median             0.0012
## Arithmetic Mean    0.0009
## Geometric Mean     0.0006
## Quartile 3         0.0125
## Maximum            0.1397
## SE Mean            0.0004
## LCL Mean (0.95)    0.0001
## UCL Mean (0.95)    0.0017
## Variance           0.0006
## Stdev              0.0235
## Skewness          -0.1070
## Kurtosis           5.9038

Plot of gun portfolio returns

Plot showing portfolio returns and the dates of mass shootings (vertical dash lines) on function of time

gun_portfolio_returns_plot <- ggplot(data = gun_co_portfolio_returns, aes(x=date, y=Ra))
gun_portfolio_returns_plot <- gun_portfolio_returns_plot + geom_line(aes(color="red"), size=1)
gun_portfolio_returns_plot <- gun_portfolio_returns_plot + theme(legend.position = "none")

## Vertical lines to show mass shootings
gun_portfolio_returns_plot <- gun_portfolio_returns_plot + geom_vline(data=us_mass03, aes(xintercept=date), linetype=3)

gun_portfolio_returns_plot <- gun_portfolio_returns_plot + labs(title = "Returns of US Gun Manufacture Portfolio", 
                                                                subtitle = "Mass Shootings with seven or more total victims", 
                                                                caption = "Source: Yahoo finance and Kaggle's Mass Shooting in United States (2018-2022)")
gun_portfolio_returns_plot <- gun_portfolio_returns_plot + ylab("Ra") + ylim(-0.2, 0.2)
gun_portfolio_returns_plot <- gun_portfolio_returns_plot + scale_x_date(date_breaks = "2 years", date_labels = "%Y")
gun_portfolio_returns_plot <- gun_portfolio_returns_plot + geom_hline(yintercept = 0, color="black")

gun_portfolio_returns_plot

ggsave("89plot09B.png")
## Saving 7 x 5 in image

Bar plot of portfolio returns

Bar plot showing the daily return of gun industry portfolio

gun_portfolio_returns_bars <- ggplot(data = gun_co_portfolio_returns, aes(x=date, y=Ra))
gun_portfolio_returns_bars <- gun_portfolio_returns_bars + geom_bar(stat = "identity", fill=palette_light()[[1]])
gun_portfolio_returns_bars <- gun_portfolio_returns_bars + geom_smooth(method = "lm")
gun_portfolio_returns_bars <- gun_portfolio_returns_bars + theme(legend.position = "none")

## Vertical lines to show mass shootings
gun_portfolio_returns_bars <- gun_portfolio_returns_bars + geom_vline(data=us_mass03, aes(xintercept=date), linetype=3)

gun_portfolio_returns_bars <- gun_portfolio_returns_bars + labs(title = "Returns of US Gun Manufacture Portfolio", 
                                                                subtitle = "Mass Shootings with seven or more total victims", 
                                                                caption = "Source: Yahoo finance and Kaggle's Mass Shooting in United States (2018-2022)")
gun_portfolio_returns_bars <- gun_portfolio_returns_bars + ylab("Ra") + ylim(-0.05,0.05)
gun_portfolio_returns_bars <- gun_portfolio_returns_bars + scale_x_date(date_breaks = "2 years", date_labels = "%Y")
gun_portfolio_returns_bars <- gun_portfolio_returns_bars + geom_hline(yintercept = 0, color="black")
gun_portfolio_returns_bars
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 140 rows containing non-finite values (stat_smooth).
## Warning: Removed 140 rows containing missing values (position_stack).

ggsave("89plot10B.png")
## Saving 7 x 5 in image
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 140 rows containing non-finite values (stat_smooth).
## Removed 140 rows containing missing values (position_stack).

plot of portfolio drawdowns

Plot showing portfolio drawdowns and mass shootings (vertical das lines) on function of time

gun_portfolio_DD_plot <- ggplot(data = gun_co_portfolio_returns, aes(x=date, y=gunDD))
gun_portfolio_DD_plot <- gun_portfolio_DD_plot + geom_line(aes(color="red"), size=1)
gun_portfolio_DD_plot <- gun_portfolio_DD_plot + theme(legend.position = "none")

## Vertical lines to show mass shootings
gun_portfolio_DD_plot <- gun_portfolio_DD_plot + geom_vline(data=us_mass03, aes(xintercept=date), linetype=3)

gun_portfolio_DD_plot <- gun_portfolio_DD_plot + labs(title = "DrawDowns of US Gun Manufacture Portfolio", 
                                                                subtitle = "Mass Shootings with seven or more total victims", 
                                                                caption = "Source: Yahoo finance and Kaggle's Mass Shooting in United States (2018-2022)")
gun_portfolio_DD_plot <- gun_portfolio_DD_plot + ylab("Draw Downs") + ylim(-0.5, 0)
gun_portfolio_DD_plot <- gun_portfolio_DD_plot + scale_x_date(date_breaks = "2 years", date_labels = "%Y")

gun_portfolio_DD_plot

ggsave("89plot11B.png")
## Saving 7 x 5 in image

Portfolio return moving averages

Plot showing 5 day moving average and mass shooting events on function of time

gun_portfolio_RaAvg <- ggplot(data = gun_co_portfolio_returns, aes(x=date, y=Ra_5day))
gun_portfolio_RaAvg <- gun_portfolio_RaAvg + geom_line(aes(color="red"), size=1)
gun_portfolio_RaAvg <- gun_portfolio_RaAvg + theme(legend.position = "none")

## Vertical lines to show mass shootings
gun_portfolio_RaAvg <- gun_portfolio_RaAvg + geom_vline(data=us_mass03, aes(xintercept=date), linetype=3)

gun_portfolio_RaAvg <- gun_portfolio_RaAvg + labs(title = "5 day moving average of US Gun Manufacture Portfolio", 
                                                                subtitle = "Mass Shootings with seven or more total victims", 
                                                                caption = "Source: Yahoo finance and Kaggle's Mass Shooting in United States (2018-2022)")
gun_portfolio_RaAvg <- gun_portfolio_RaAvg + ylab("Draw Downs")
gun_portfolio_RaAvg <- gun_portfolio_RaAvg + scale_x_date(date_breaks = "2 years", date_labels = "%Y")

gun_portfolio_RaAvg

ggsave("89plot12B.png")
## Saving 7 x 5 in image

Portfolio return momentum

Plot showing 10 day portfolio momentun and mass shootign events on function of time

gun_portfolio_RaAvg <- ggplot(data = gun_co_portfolio_returns, aes(x=date, y=Mom_10day))
gun_portfolio_RaAvg <- gun_portfolio_RaAvg + geom_line(aes(color="red"), size=1)
gun_portfolio_RaAvg <- gun_portfolio_RaAvg + theme(legend.position = "none")

## Vertical lines to show mass shootings
gun_portfolio_RaAvg <- gun_portfolio_RaAvg + geom_vline(data=us_mass03, aes(xintercept=date), linetype=3)

gun_portfolio_RaAvg <- gun_portfolio_RaAvg + labs(title = "10 day return momentum of US Gun Manufacture Portfolio", 
                                                                subtitle = "Mass Shootings with seven or more total victims", 
                                                                caption = "Source: Yahoo finance and Kaggle's Mass Shooting in United States (2018-2022)")
gun_portfolio_RaAvg <- gun_portfolio_RaAvg + ylab("Draw Downs")
gun_portfolio_RaAvg <- gun_portfolio_RaAvg + scale_x_date(date_breaks = "2 years", date_labels = "%Y")

gun_portfolio_RaAvg
## Warning: Removed 10 row(s) containing missing values (geom_path).

ggsave("89plot13B.png")
## Saving 7 x 5 in image
## Warning: Removed 10 row(s) containing missing values (geom_path).

Portfolio’s Cumulative return

Cumulative return using tq_portfolio function

gun_co_portfolio_growth <- gun_co_returns %>% tq_portfolio(assets_col = symbol, returns_col = Ra, weights = gun_wts, col_rename = "investment.growth", wealth.index=TRUE)
gun_co_portfolio_growth$investment.growth <- gun_co_portfolio_growth$investment.growth-1
head(gun_co_portfolio_growth)
## # A tibble: 6 × 2
##   date       investment.growth
##   <date>                 <dbl>
## 1 2009-10-01        0         
## 2 2009-10-02       -0.0272    
## 3 2009-10-05       -0.0350    
## 4 2009-10-06       -0.0163    
## 5 2009-10-07       -0.00777   
## 6 2009-10-08       -0.00000128

Portfolio’s compounded return plot

Plot showing cumulative portfolio return and mass shooitng events on function of time

gun_portfolio_cum_return <- ggplot(data = gun_co_portfolio_growth, aes(x=date, y=investment.growth))
gun_portfolio_cum_return <- gun_portfolio_cum_return + geom_line(size=1, color=palette_light()[[1]])
gun_portfolio_cum_return <- gun_portfolio_cum_return + geom_smooth(method = "loess")
gun_portfolio_cum_return <- gun_portfolio_cum_return + theme(legend.position = "none")

## Vertical lines to show mass shootings
gun_portfolio_cum_return <- gun_portfolio_cum_return + geom_vline(data=us_mass03, aes(xintercept=date), linetype=3)

gun_portfolio_cum_return <- gun_portfolio_cum_return + labs(title = "Compounded Return of US Gun Manufacture Portfolio", 
                                                            subtitle = "Mass Shootings with seven or more total victims", 
                                                            caption = "Source: Yahoo finance and Kaggle's Mass Shooting in United States (2018-2022)")
gun_portfolio_cum_return <- gun_portfolio_cum_return + ylab("Ra") 
gun_portfolio_cum_return <- gun_portfolio_cum_return + scale_x_date(date_breaks = "2 years", date_labels = "%Y")

gun_portfolio_cum_return
## `geom_smooth()` using formula 'y ~ x'

ggsave("89plot14B.png")
## Saving 7 x 5 in image
## `geom_smooth()` using formula 'y ~ x'

Join portfolio dataframe with mass shooting dataframe

Merging mass shooting data set and portfolio data set. Those observations where 20 day moving average is NA were eliminated. An indicator (dummy) variable for shooting was created

temp_mass_dates <- us_mass02 %>% filter(total_victims >= 7) %>% select("date", "fatalities", "injured", "total_victims")
merged_df <- merge(gun_co_portfolio_returns, temp_mass_dates, by.x = "date", by.y = "date", all.x = TRUE)
merged_df <- merged_df %>%  replace_na(list(fatalities=0, injured=0, total_victims=0)) %>%
                        mutate(shooting = ifelse(total_victims>0,1,0))
merged_df2 <- merged_df %>% drop_na(Mom_20day) #%>% drop_na(Ra_10day)
tail(merged_df2, 10)
##            date            Ra      SMA.2day    SMA.15day      gunDD gunDD_2day
## 3132 2022-06-09  4.580987e-03 -0.0041478253  0.004347252 -0.3335843 -0.3554429
## 3133 2022-06-10 -2.838113e-03  0.0008714373  0.004435712 -0.3354756 -0.3554188
## 3134 2022-06-13 -3.004751e-02 -0.0164428100  0.001333129 -0.3554429 -0.3522838
## 3135 2022-06-14  3.746473e-05 -0.0150050213  0.001744966 -0.3554188 -0.3804124
## 3136 2022-06-15  4.863690e-03  0.0024505774  0.001966289 -0.3522838 -0.3838896
## 3137 2022-06-16 -4.342748e-02 -0.0192818958 -0.004065449 -0.3804124 -0.3841058
## 3138 2022-06-17 -5.611991e-03 -0.0245197364 -0.005448874 -0.3838896 -0.3927233
## 3139 2022-06-21 -3.510373e-04 -0.0029815143 -0.005471049 -0.3841058 -0.3599289
## 3140 2022-06-22 -1.399185e-02 -0.0071714430 -0.007230260 -0.3927233 -0.3262850
## 3141 2022-06-23  5.400243e-02  0.0200052899 -0.003036847 -0.3599289 -0.3580488
##      gunDD_3day gunDD_5day gunDD_10day gunDD_15day gunDD_20day       Ra_2day
## 3132 -0.3554188 -0.3804124  -0.3262850  -0.3679977  -0.3731047 -0.0041478253
## 3133 -0.3522838 -0.3838896  -0.3580488  -0.3491225  -0.3785227  0.0008714373
## 3134 -0.3804124 -0.3841058  -0.3635919  -0.3760208  -0.3752875 -0.0164428100
## 3135 -0.3838896 -0.3927233  -0.3728398  -0.3685268  -0.3863849 -0.0150050213
## 3136 -0.3841058 -0.3599289  -0.3768972  -0.3757678  -0.3820448  0.0024505774
## 3137 -0.3927233 -0.3262850  -0.3679977  -0.3731047  -0.3950705 -0.0192818958
## 3138 -0.3599289 -0.3580488  -0.3491225  -0.3785227  -0.3872490 -0.0245197364
## 3139 -0.3262850 -0.3635919  -0.3760208  -0.3752875  -0.3773754 -0.0029815143
## 3140 -0.3580488 -0.3728398  -0.3685268  -0.3863849  -0.3727382 -0.0071714430
## 3141 -0.3635919 -0.3768972  -0.3757678  -0.3820448  -0.3760701  0.0200052899
##            Ra_3day      Ra_5day     Ra_10day     Mom_2day     Mom_3day
## 3132 -4.185704e-05 -0.008620005  0.001043022 -0.034628495 -0.004543523
## 3133 -3.711254e-03 -0.003272175 -0.000754728  0.002875577  0.007701803
## 3134 -9.434878e-03 -0.006602238 -0.003757638  0.034911197 -0.013379974
## 3135 -1.094939e-02 -0.008228761 -0.004993524 -0.043464946 -0.005649456
## 3136 -8.382117e-03 -0.004680696 -0.003617277 -0.010475681 -0.005214727
## 3137 -1.284211e-02 -0.014282389 -0.011451197  0.043076444  0.029435633
## 3138 -1.472526e-02 -0.014837165 -0.009054670 -0.008379857  0.059614420
## 3139 -1.646350e-02 -0.008897871 -0.007750055  0.054353466  0.052913839
## 3140 -6.651626e-03 -0.011703734 -0.009966247  0.066554650 -0.033155442
## 3141  1.321985e-02 -0.001875986 -0.003278341 -0.101149720 -0.062637205
##           Mom_5day   Mom_10day    Mom_15day     Mom_20day fatalities injured
## 3132 -0.0480084690  0.04798181  0.009701543 -0.0003147059          0       0
## 3133 -0.0027738785 -0.04430918  0.032703896 -0.0058045291          0       0
## 3134  0.0296964700  0.02141273 -0.011278637  0.0352531855          0       0
## 3135 -0.0140293134 -0.01456875  0.011972425 -0.0178014742          0       0
## 3136  0.0491387385 -0.01133326 -0.016330538  0.0022093910          0       0
## 3137  0.0959902830  0.05771001  0.047693763  0.0223487799          0       0
## 3138 -0.0415352998  0.03547777 -0.003030651  0.0185414733          0       0
## 3139 -0.0082837389 -0.04097511  0.005556716  0.0164646274          0       0
## 3140 -0.0005394359  0.02600174 -0.003772161  0.0214397378          0       0
## 3141 -0.0604720002 -0.06546928 -0.046929348 -0.0593143249          0       0
##      total_victims shooting
## 3132             0        0
## 3133             0        0
## 3134             0        0
## 3135             0        0
## 3136             0        0
## 3137             0        0
## 3138             0        0
## 3139             0        0
## 3140             0        0
## 3141             0        0

correlation plot

A correlation plot was set up in order to determine multicollinearity among variables

corrplot::corrplot(cor(merged_df2[,-1]),  method="number",  number.cex = 0.7) #addCoef.col=1,

ggsave("89plot15B.png")
## Saving 7 x 5 in image
## `geom_smooth()` using formula 'y ~ x'

Splitting data set in trainig and test set

The data set was split into train and test. Until December 31, 2019 was used as a training, after that was used for testing

merged_df_train <- merged_df2 %>% filter(date<="2019-12-31")

merged_df_test <- merged_df2 %>% filter(date > "2019-12-31")
tail(merged_df_test, 10)
##           date            Ra      SMA.2day    SMA.15day      gunDD gunDD_2day
## 615 2022-06-09  4.580987e-03 -0.0041478253  0.004347252 -0.3335843 -0.3554429
## 616 2022-06-10 -2.838113e-03  0.0008714373  0.004435712 -0.3354756 -0.3554188
## 617 2022-06-13 -3.004751e-02 -0.0164428100  0.001333129 -0.3554429 -0.3522838
## 618 2022-06-14  3.746473e-05 -0.0150050213  0.001744966 -0.3554188 -0.3804124
## 619 2022-06-15  4.863690e-03  0.0024505774  0.001966289 -0.3522838 -0.3838896
## 620 2022-06-16 -4.342748e-02 -0.0192818958 -0.004065449 -0.3804124 -0.3841058
## 621 2022-06-17 -5.611991e-03 -0.0245197364 -0.005448874 -0.3838896 -0.3927233
## 622 2022-06-21 -3.510373e-04 -0.0029815143 -0.005471049 -0.3841058 -0.3599289
## 623 2022-06-22 -1.399185e-02 -0.0071714430 -0.007230260 -0.3927233 -0.3262850
## 624 2022-06-23  5.400243e-02  0.0200052899 -0.003036847 -0.3599289 -0.3580488
##     gunDD_3day gunDD_5day gunDD_10day gunDD_15day gunDD_20day       Ra_2day
## 615 -0.3554188 -0.3804124  -0.3262850  -0.3679977  -0.3731047 -0.0041478253
## 616 -0.3522838 -0.3838896  -0.3580488  -0.3491225  -0.3785227  0.0008714373
## 617 -0.3804124 -0.3841058  -0.3635919  -0.3760208  -0.3752875 -0.0164428100
## 618 -0.3838896 -0.3927233  -0.3728398  -0.3685268  -0.3863849 -0.0150050213
## 619 -0.3841058 -0.3599289  -0.3768972  -0.3757678  -0.3820448  0.0024505774
## 620 -0.3927233 -0.3262850  -0.3679977  -0.3731047  -0.3950705 -0.0192818958
## 621 -0.3599289 -0.3580488  -0.3491225  -0.3785227  -0.3872490 -0.0245197364
## 622 -0.3262850 -0.3635919  -0.3760208  -0.3752875  -0.3773754 -0.0029815143
## 623 -0.3580488 -0.3728398  -0.3685268  -0.3863849  -0.3727382 -0.0071714430
## 624 -0.3635919 -0.3768972  -0.3757678  -0.3820448  -0.3760701  0.0200052899
##           Ra_3day      Ra_5day     Ra_10day     Mom_2day     Mom_3day
## 615 -4.185704e-05 -0.008620005  0.001043022 -0.034628495 -0.004543523
## 616 -3.711254e-03 -0.003272175 -0.000754728  0.002875577  0.007701803
## 617 -9.434878e-03 -0.006602238 -0.003757638  0.034911197 -0.013379974
## 618 -1.094939e-02 -0.008228761 -0.004993524 -0.043464946 -0.005649456
## 619 -8.382117e-03 -0.004680696 -0.003617277 -0.010475681 -0.005214727
## 620 -1.284211e-02 -0.014282389 -0.011451197  0.043076444  0.029435633
## 621 -1.472526e-02 -0.014837165 -0.009054670 -0.008379857  0.059614420
## 622 -1.646350e-02 -0.008897871 -0.007750055  0.054353466  0.052913839
## 623 -6.651626e-03 -0.011703734 -0.009966247  0.066554650 -0.033155442
## 624  1.321985e-02 -0.001875986 -0.003278341 -0.101149720 -0.062637205
##          Mom_5day   Mom_10day    Mom_15day     Mom_20day fatalities injured
## 615 -0.0480084690  0.04798181  0.009701543 -0.0003147059          0       0
## 616 -0.0027738785 -0.04430918  0.032703896 -0.0058045291          0       0
## 617  0.0296964700  0.02141273 -0.011278637  0.0352531855          0       0
## 618 -0.0140293134 -0.01456875  0.011972425 -0.0178014742          0       0
## 619  0.0491387385 -0.01133326 -0.016330538  0.0022093910          0       0
## 620  0.0959902830  0.05771001  0.047693763  0.0223487799          0       0
## 621 -0.0415352998  0.03547777 -0.003030651  0.0185414733          0       0
## 622 -0.0082837389 -0.04097511  0.005556716  0.0164646274          0       0
## 623 -0.0005394359  0.02600174 -0.003772161  0.0214397378          0       0
## 624 -0.0604720002 -0.06546928 -0.046929348 -0.0593143249          0       0
##     total_victims shooting
## 615             0        0
## 616             0        0
## 617             0        0
## 618             0        0
## 619             0        0
## 620             0        0
## 621             0        0
## 622             0        0
## 623             0        0
## 624             0        0

Linear regressions

The following lines are the series of models built for this analysis. Transformations were implemented as well

SMA Short Model 1
Model1 <- lm(data = merged_df_train, formula = SMA.2day ~ fatalities + injured)
summary(Model1)
## 
## Call:
## lm(formula = SMA.2day ~ fatalities + injured, data = merged_df_train)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.129373 -0.008098  0.000027  0.008098  0.086079 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)   
## (Intercept) 8.606e-04  3.226e-04   2.667  0.00769 **
## fatalities  3.220e-05  2.240e-04   0.144  0.88570   
## injured     1.756e-05  4.001e-05   0.439  0.66083   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01611 on 2514 degrees of freedom
## Multiple R-squared:  0.0002286,  Adjusted R-squared:  -0.0005668 
## F-statistic: 0.2874 on 2 and 2514 DF,  p-value: 0.7502
plot(Model1)

## Warning in sqrt(crit * p * (1 - hh)/hh): NaNs produced

## Warning in sqrt(crit * p * (1 - hh)/hh): NaNs produced

Model 5
Model5 <- lm(data = merged_df_train, formula = SMA.2day ~ total_victims)
summary(Model5)
## 
## Call:
## lm(formula = SMA.2day ~ total_victims, data = merged_df_train)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.129375 -0.008100  0.000025  0.008096  0.086077 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   8.621e-04  3.214e-04   2.683  0.00735 **
## total_victims 1.934e-05  2.557e-05   0.756  0.44967   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01611 on 2515 degrees of freedom
## Multiple R-squared:  0.0002273,  Adjusted R-squared:  -0.0001703 
## F-statistic: 0.5717 on 1 and 2515 DF,  p-value: 0.4497
plot(Model5)

SMA long vs tatalities + injured Model 2
Model2 <- lm(data = merged_df_train, formula = SMA.15day ~ fatalities + injured)
summary(Model2)
## 
## Call:
## lm(formula = SMA.15day ~ fatalities + injured, data = merged_df_train)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0219520 -0.0039910  0.0001739  0.0037688  0.0198369 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  8.870e-04  1.165e-04   7.614 3.73e-14 ***
## fatalities  -1.879e-04  8.089e-05  -2.323  0.02027 *  
## injured      4.036e-05  1.445e-05   2.793  0.00526 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.005818 on 2514 degrees of freedom
## Multiple R-squared:  0.003211,   Adjusted R-squared:  0.002418 
## F-statistic: 4.049 on 2 and 2514 DF,  p-value: 0.01756
par(mfrow=c(1,2))
plot(Model2)

## Warning in sqrt(crit * p * (1 - hh)/hh): NaNs produced

## Warning in sqrt(crit * p * (1 - hh)/hh): NaNs produced

require(ggiraph)
## Loading required package: ggiraph
## Warning: package 'ggiraph' was built under R version 4.2.1
require(ggiraphExtra)
## Loading required package: ggiraphExtra
## Warning: package 'ggiraphExtra' was built under R version 4.2.1
require(plyr)
## Loading required package: plyr
## ------------------------------------------------------------------------------
## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)
## ------------------------------------------------------------------------------
## 
## Attaching package: 'plyr'
## The following object is masked from 'package:purrr':
## 
##     compact
## The following objects are masked from 'package:dplyr':
## 
##     arrange, count, desc, failwith, id, mutate, rename, summarise,
##     summarize
#install.packages("devtools")
#devtools::install_github("cardiomoon/ggiraphExtra")
#install.packages("predict3d")
ggPredict(Model2,se=TRUE,interactive=TRUE)
SMA long vs total victims Model 6
Model6 <- lm(data = merged_df_train, formula = SMA.15day ~ total_victims)
summary(Model6)
## 
## Call:
## lm(formula = SMA.15day ~ total_victims, data = merged_df_train)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0219273 -0.0039855  0.0001545  0.0037848  0.0198616 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   8.624e-04  1.162e-04   7.421 1.58e-13 ***
## total_victims 1.264e-05  9.247e-06   1.366    0.172    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.005824 on 2515 degrees of freedom
## Multiple R-squared:  0.0007419,  Adjusted R-squared:  0.0003446 
## F-statistic: 1.867 on 1 and 2515 DF,  p-value: 0.1719
plot(Model6)

Draw Down 2 day
DD1 <- lm(data = merged_df_train, formula = gunDD_2day ~ total_victims)
summary(DD1)
## 
## Call:
## lm(formula = gunDD_2day ~ total_victims, data = merged_df_train)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.33012 -0.11365  0.01923  0.10485  0.22278 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -0.2175150  0.0027745 -78.397   <2e-16 ***
## total_victims -0.0002633  0.0002208  -1.193    0.233    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1391 on 2515 degrees of freedom
## Multiple R-squared:  0.0005654,  Adjusted R-squared:  0.000168 
## F-statistic: 1.423 on 1 and 2515 DF,  p-value: 0.2331
Draw Down 3 day
DD2 <- lm(data = merged_df_train, formula = gunDD_3day ~ total_victims)
summary(DD2)
## 
## Call:
## lm(formula = gunDD_3day ~ total_victims, data = merged_df_train)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.33002 -0.11396  0.01924  0.10494  0.22180 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -0.2176200  0.0027756 -78.404   <2e-16 ***
## total_victims -0.0002459  0.0002209  -1.113    0.266    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1391 on 2515 degrees of freedom
## Multiple R-squared:  0.0004924,  Adjusted R-squared:  9.498e-05 
## F-statistic: 1.239 on 1 and 2515 DF,  p-value: 0.2658
Draw Down 5 day
DD3 <- lm(data = merged_df_train, formula = gunDD_5day ~ total_victims)
summary(DD3)
## 
## Call:
## lm(formula = gunDD_5day ~ total_victims, data = merged_df_train)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.32985 -0.11403  0.01925  0.10523  0.22367 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -0.2177914  0.0027778 -78.405   <2e-16 ***
## total_victims -0.0002937  0.0002210  -1.329    0.184    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1392 on 2515 degrees of freedom
## Multiple R-squared:  0.0007016,  Adjusted R-squared:  0.0003043 
## F-statistic: 1.766 on 1 and 2515 DF,  p-value: 0.184
Draw Down 10 day
DD4 <- lm(data = merged_df_train, formula = gunDD_10day ~ total_victims)
summary(DD4)
## 
## Call:
## lm(formula = gunDD_10day ~ total_victims, data = merged_df_train)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.32939 -0.11425  0.01935  0.10557  0.21825 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -0.2182529  0.0027830 -78.425   <2e-16 ***
## total_victims -0.0003128  0.0002215  -1.413    0.158    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1395 on 2515 degrees of freedom
## Multiple R-squared:  0.0007928,  Adjusted R-squared:  0.0003955 
## F-statistic: 1.995 on 1 and 2515 DF,  p-value: 0.1579
Draw Down 15 day
DD5 <- lm(data = merged_df_train, formula = gunDD_15day ~ total_victims)
summary(DD5)
## 
## Call:
## lm(formula = gunDD_15day ~ total_victims, data = merged_df_train)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.32888 -0.11420  0.01863  0.10624  0.22320 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -0.2187592  0.0027882 -78.459   <2e-16 ***
## total_victims -0.0002610  0.0002219  -1.176     0.24    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1397 on 2515 degrees of freedom
## Multiple R-squared:  0.0005498,  Adjusted R-squared:  0.0001524 
## F-statistic: 1.383 on 1 and 2515 DF,  p-value: 0.2396
Draw Down 20 day
DD6 <- lm(data = merged_df_train, formula = gunDD_20day ~ total_victims)
summary(DD6)
## 
## Call:
## lm(formula = gunDD_20day ~ total_victims, data = merged_df_train)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.32846 -0.11512  0.01825  0.10651  0.21918 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -0.2191758  0.0027932 -78.467   <2e-16 ***
## total_victims -0.0002939  0.0002223  -1.322    0.186    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.14 on 2515 degrees of freedom
## Multiple R-squared:  0.0006949,  Adjusted R-squared:  0.0002976 
## F-statistic: 1.749 on 1 and 2515 DF,  p-value: 0.1861
Momentum 2 day
mom1 <- lm(data = merged_df_train, formula = Mom_2day ~ total_victims)
summary(mom1)
## 
## Call:
## lm(formula = Mom_2day ~ total_victims, data = merged_df_train)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.19071 -0.01610 -0.00056  0.01634  0.19255 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)  
## (Intercept)    3.598e-05  6.384e-04   0.056   0.9551  
## total_victims -1.147e-04  5.080e-05  -2.257   0.0241 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.03199 on 2515 degrees of freedom
## Multiple R-squared:  0.002022,   Adjusted R-squared:  0.001625 
## F-statistic: 5.095 on 1 and 2515 DF,  p-value: 0.02409
Momentum 3 day
mom2 <- lm(data = merged_df_train, formula = Mom_3day ~ total_victims)
summary(mom2)
## 
## Call:
## lm(formula = Mom_3day ~ total_victims, data = merged_df_train)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.170502 -0.017780 -0.000561  0.017318  0.229605 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)
## (Intercept)    1.330e-05  6.551e-04   0.020    0.984
## total_victims -4.235e-05  5.213e-05  -0.812    0.417
## 
## Residual standard error: 0.03283 on 2515 degrees of freedom
## Multiple R-squared:  0.0002624,  Adjusted R-squared:  -0.0001351 
## F-statistic:  0.66 on 1 and 2515 DF,  p-value: 0.4166
Momentum 5 day vs fatalities + injured Model 3
Model3 <- lm(data = merged_df_train, formula = Mom_5day ~ fatalities + injured)
summary(Model3)
## 
## Call:
## lm(formula = Mom_5day ~ fatalities + injured, data = merged_df_train)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.169070 -0.016474 -0.000299  0.017039  0.191508 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept)  5.270e-05  6.530e-04   0.081    0.936
## fatalities  -4.435e-05  4.534e-04  -0.098    0.922
## injured     -1.295e-04  8.099e-05  -1.599    0.110
## 
## Residual standard error: 0.03261 on 2514 degrees of freedom
## Multiple R-squared:  0.002112,   Adjusted R-squared:  0.001318 
## F-statistic: 2.661 on 2 and 2514 DF,  p-value: 0.0701
plot(Model3)

## Warning in sqrt(crit * p * (1 - hh)/hh): NaNs produced

## Warning in sqrt(crit * p * (1 - hh)/hh): NaNs produced

Momentum 5 day vs total victims Model 7
Model7 <- lm(data = merged_df_train, formula = Mom_5day ~ total_victims)
summary(Model7)
## 
## Call:
## lm(formula = Mom_5day ~ total_victims, data = merged_df_train)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.169079 -0.016358 -0.000308  0.017030  0.191499 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)  
## (Intercept)    6.189e-05  6.505e-04   0.095   0.9242  
## total_victims -1.191e-04  5.177e-05  -2.301   0.0215 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0326 on 2515 degrees of freedom
## Multiple R-squared:  0.002101,   Adjusted R-squared:  0.001704 
## F-statistic: 5.296 on 1 and 2515 DF,  p-value: 0.02146
par(mfrow=c(2,2))
plot(Model7)

ggPredict(Model7,se=TRUE,interactive=TRUE)
Momentum 10 day
mom4 <- lm(data = merged_df_train, formula = Mom_10day ~ total_victims)
summary(mom4)
## 
## Call:
## lm(formula = Mom_10day ~ total_victims, data = merged_df_train)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.174663 -0.017609  0.000211  0.017679  0.163871 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)
## (Intercept)    2.094e-05  6.423e-04   0.033    0.974
## total_victims -4.767e-05  5.111e-05  -0.933    0.351
## 
## Residual standard error: 0.03219 on 2515 degrees of freedom
## Multiple R-squared:  0.0003457,  Adjusted R-squared:  -5.176e-05 
## F-statistic: 0.8698 on 1 and 2515 DF,  p-value: 0.3511
Momentum 15 day
mom5 <- lm(data = merged_df_train, formula = Mom_15day ~ total_victims )
summary(mom5)
## 
## Call:
## lm(formula = Mom_15day ~ total_victims, data = merged_df_train)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.171799 -0.017124  0.000485  0.017619  0.191330 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)  
## (Intercept)    4.578e-05  6.491e-04   0.071   0.9438  
## total_victims -9.134e-05  5.165e-05  -1.768   0.0771 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.03253 on 2515 degrees of freedom
## Multiple R-squared:  0.001242,   Adjusted R-squared:  0.0008448 
## F-statistic: 3.127 on 1 and 2515 DF,  p-value: 0.07711
Momentum 20 day vs Fatalities + Injured linear - linear Model 4
Model4 <- lm(data = merged_df_train, formula = Mom_20day ~ fatalities + injured)
summary(Model4)
## 
## Call:
## lm(formula = Mom_20day ~ fatalities + injured, data = merged_df_train)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.181020 -0.017133 -0.000324  0.016720  0.206396 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept)  1.075e-04  6.569e-04   0.164    0.870
## fatalities  -4.548e-04  4.561e-04  -0.997    0.319
## injured     -6.589e-05  8.147e-05  -0.809    0.419
## 
## Residual standard error: 0.0328 on 2514 degrees of freedom
## Multiple R-squared:  0.002098,   Adjusted R-squared:  0.001304 
## F-statistic: 2.643 on 2 and 2514 DF,  p-value: 0.07134
plot(Model4)

## Warning in sqrt(crit * p * (1 - hh)/hh): NaNs produced

## Warning in sqrt(crit * p * (1 - hh)/hh): NaNs produced

Momentum 20 day vs Total Victims linear - linear Model 8
Model8 <- lm(data = merged_df_train, formula = Mom_20day ~ total_victims)
summary(Model8)
## 
## Call:
## lm(formula = Mom_20day ~ total_victims, data = merged_df_train)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.180978 -0.017098 -0.000289  0.016761  0.206438 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)  
## (Intercept)    6.552e-05  6.545e-04   0.100   0.9203  
## total_victims -1.131e-04  5.208e-05  -2.172   0.0299 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0328 on 2515 degrees of freedom
## Multiple R-squared:  0.001873,   Adjusted R-squared:  0.001476 
## F-statistic: 4.718 on 1 and 2515 DF,  p-value: 0.02993
plot(Model8)

Momentum 20 day linear - log Model 8a
Model8a <- lm(data = merged_df_train, formula = Mom_20day ~ log(total_victims+1))
summary(Model8a)
## 
## Call:
## lm(formula = Mom_20day ~ log(total_victims + 1), data = merged_df_train)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.181023 -0.017137 -0.000328  0.016766  0.206392 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)
## (Intercept)             0.0001111  0.0006592   0.169    0.866
## log(total_victims + 1) -0.0022828  0.0016907  -1.350    0.177
## 
## Residual standard error: 0.03282 on 2515 degrees of freedom
## Multiple R-squared:  0.0007243,  Adjusted R-squared:  0.000327 
## F-statistic: 1.823 on 1 and 2515 DF,  p-value: 0.1771
Momentum 20 day log - linear Model 8b
c <- abs(min(merged_df_train$Mom_20day))
c
## [1] 0.180912
merged_df_train$log_mom_20day_c <- log(merged_df_train$Mom_20day+c+1)

Model8b <- lm(data = merged_df_train, formula = log_mom_20day_c ~ total_victims)
summary(Model8b)
## 
## Call:
## lm(formula = log_mom_20day_c ~ total_victims, data = merged_df_train)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.165957 -0.014198  0.000141  0.014478  0.161486 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    1.660e-01  5.556e-04 298.716   <2e-16 ***
## total_victims -9.839e-05  4.421e-05  -2.225   0.0261 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.02784 on 2515 degrees of freedom
## Multiple R-squared:  0.001965,   Adjusted R-squared:  0.001569 
## F-statistic: 4.953 on 1 and 2515 DF,  p-value: 0.02614
par(mfrow=c(1,2))
plot(Model8b)

ggPredict(Model8b, se=TRUE,interactive=TRUE)
Momentum 20 day log - log Model 8c
c <- abs(min(merged_df_train$Mom_20day))
c
## [1] 0.180912
Model8c <- lm(data = merged_df_train, formula = log(Mom_20day+c+1) ~ log(total_victims+1))
summary(Model8c)
## 
## Call:
## lm(formula = log(Mom_20day + c + 1) ~ log(total_victims + 1), 
##     data = merged_df_train)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.166000 -0.014235  0.000104  0.014478  0.161443 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             0.1659997  0.0005596  296.62   <2e-16 ***
## log(total_victims + 1) -0.0020517  0.0014352   -1.43    0.153    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.02786 on 2515 degrees of freedom
## Multiple R-squared:  0.0008119,  Adjusted R-squared:  0.0004146 
## F-statistic: 2.044 on 1 and 2515 DF,  p-value: 0.153

Prediction - Model 2

library(forecast)
predict_01 <- predict(Model2, newdata = merged_df_test)
head(predict_01,10)
##            1            2            3            4            5            6 
## 0.0008870347 0.0008870347 0.0008870347 0.0008870347 0.0008870347 0.0008870347 
##            7            8            9           10 
## 0.0008870347 0.0008870347 0.0008870347 0.0008870347
accuracy_01 <- accuracy(predict_01, merged_df_test$SMA.15day)
accuracy_01
##                   ME        RMSE         MAE      MPE     MAPE
## Test set 0.000321654 0.005971153 0.004634279 82.64277 140.5359

Prediction - Model 7

predict_02 <- predict(Model7, newdata = merged_df_test)
head(predict_02,10)
##            1            2            3            4            5            6 
## 6.189288e-05 6.189288e-05 6.189288e-05 6.189288e-05 6.189288e-05 6.189288e-05 
##            7            8            9           10 
## 6.189288e-05 6.189288e-05 6.189288e-05 6.189288e-05
accuracy_02 <- accuracy(predict_02, merged_df_test$Mom_5day)
accuracy_02
##                     ME       RMSE        MAE      MPE     MAPE
## Test set -0.0001480064 0.03681865 0.02652235 93.54564 106.0142

The next two plots were generated so we can confirm that our code was corrected

Calculating compound return of gun portfolio

gun_co_portfolio_returns1 <- gun_co_portfolio_returns
gun_co_portfolio_returns1 <- xts(gun_co_portfolio_returns1[,-1], order.by = gun_co_portfolio_returns$date,)
#head(gun_co_portfolio_returns1)

gun_comp_return <- Return.cumulative(gun_co_portfolio_returns1$Ra, geometric = TRUE)
gun_comp_return
##                         Ra
## Cumulative Return 6.729115
chart.CumReturns(gun_co_portfolio_returns1$Ra, wealth.index = FALSE, geometric = TRUE)

Estimating gun Portfolio Drawdown

library(data.table)
gun_DD_plot <- chart.Drawdown(gun_co_portfolio_returns1$Ra)
gun_DD_plot

gun_DDtable <- table.Drawdowns(gun_co_portfolio_returns1$Ra, top = 7, digits = 4)
gun_DDtable
##         From     Trough         To   Depth Length To Trough Recovery
## 1 2014-01-17 2014-12-16 2016-02-26 -0.5476    531       231      300
## 2 2016-03-23 2019-09-03 2020-07-02 -0.5389   1078       868      210
## 3 2021-07-01 2022-07-18       <NA> -0.3951    268       263       NA
## 4 2012-05-03 2012-06-07 2012-11-23 -0.3630    141        25      116
## 5 2012-11-30 2012-12-18 2013-09-11 -0.3122    196        13      183
## 6 2011-08-31 2011-10-03 2011-12-21 -0.3036     79        23       56
## 7 2020-08-11 2020-11-24 2021-06-01 -0.2743    203        75      128